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Abstract 

To explain the pulsed emission of the rotation powered pulsars from radio to gamma-ray, the polar cap 
models, the slot gap models, and the outer gap models are proposed. The recent observations suggest that 
these models arc likely to co-exist in the same magnctosphcre. If so, their mutual relation is known to be 
troublesome (Harding 2009) due to the boundary conditions and the direction of the current which are 
properly assumed in each acceleration models. 

We performed a particle simulation for the global magnetospheric structure. Based on the simulation, 
we present a new picture of the global structure of the pulsar magnetosphere. It is found that a new dead 
zone is formed along the current neutral line which separates the oppositely directed current. We shall call 
this the current- neutral zone. We suggest that the polar cap accelerators and the slot gaps locate above 
the current-neutral zone, and the outer gap exist between the current neutral zone and the traditional 
dead zone. We also give an estimate of the super-rotation region. 

Key words: magnetic fields - plasmas - pulsars: general. 



1. Introduction 

The mechanism of the pulsed emission from pulsars has 
been an important unsolved problem for more than forty 
years. The recent observations of the Fermi Gamma-ray 
Space Telescope has promoted an understanding of the 
emission from gamma-ray pulsars. The pulse shapes show 
substantial diversity, but roughly 75% of the gamma-ray 
profiles have two peaks, separated by > 0.2 of rotation 
phase (Abdo et al. 2010b). Additionally, most of the 
radio-loud gamma-ray pulsars have the first gamma-ray 
peak lagging the main radio component. This indicates 
that radio and gamma-ray beams are coming from differ- 
ent emitting regions. Also, the high-confidence gamma- 
ray from a millisecond pulsar is detected first by Fermi. 
After that, it was found that most of the millisecond pul- 
sars have the features of the gamma-ray light curves and 
spectra like normal pulsars. 

Polar cap (PC) (eg, Daugherty & Harding 1996), slot 
gap (SG) (eg, Muslimov & Harding 2003; Muslimov & 
Harding 2004) and outer gap (OG) (eg, Takata ct al. 
2004; Hirotani 2008) models are promising particle accel- 
eration models to explain pulsed emission from pulsars. 
In PC models, particle acceleration and resultant gamma- 
ray emission takes place in the open field line region above 
the magnetic poles. The gamma-rays undergo the mag- 
netic pair production, which leads to super-exponential 
spectral cutoff. Then, the produced pairs screen out the 
accelerating electric field above the so-called pair forma- 
tion front (PFF) . It is generally believed that the observed 
radio emission is basically coming from upper part of the 
PC region. The SG, which is the rim of the PC, extends 



along the last-open-field line. There is no sufficient accel- 
eration to emit gamma-rays producing electron-positron 
pairs at the low altitude because the field-aligned electric 
field near the last open field line is weak. Therefore, par- 
ticle acceleration lasts up to high altitudes of many stel- 
lar radii. For this reason, the SG can provide extended 
high-energy emission above the last-open-field lines. The 
OG is located between a neighborhood of the null charge 
surface and the light cylinder along the open field lines. 
The accelerating electric field exists because of charge de- 
pletion, and photon-photon pair production occurs in the 
gap. The OG can also provide extended high energy emis- 
sion. In addition, SG and OG models have roughly simple 
exponential spectral cutoff, which is more gradual than 
that of PC models. 

SG and/or OG models can explain the double-peaked 
gamma-ray light curves (eg, Romani <fc Watters 2010; 
Venter ct al. 2009). On the other hand, PC models 
have much more difficulty reproducing the observed wide 
pulse profile because the beam size from the PC is too 
small. Furthermore, the possibility of gamma-ray emis- 
sion from the PC, which would exhibit super-exponential 
cutoff due to magnetic pair attenuation, was eliminated by 
the spectral observation such as MAGIC and Fermi (eg, 
Aliu et al. 2008; Abdo et al. 2010b). Also, there are not 
only the study of light curves and averaged spectrum but 
also that of phase-resolved spectra (eg, Takata & Chang 
2007; Harding et al. 2008). 

In PC, SG and OG models, acceleration is caused by 
an electric field parallel to the magnetic field line, and 
the accelerating electric field is calculated by the Poisson 
equation for the non-corotational potential: 
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-V 2 $ = 47r(p-p GJ ), 
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where E = — (ft x r) x B/c — V$; p and poj = ^ • 
B/(2irc)[l — (ft x r/c) 2 ] arc the charge density and 
Goldreich- Julian charge density (Goldreich & Julian 
1969), respectively; ft is the angular momentum; and r is 
the position vector. Thus, the field-aligned electric field 
En = — V<£> • B, with B unit vector along the magnetic 
field line, is formed by the difference between local charge 
density and Goldreich- Julian density. 

The boundary conditions of these models, which are re- 
quired to solve the Poisson equation, are summarized as 
follows. Most of the PC and SG models assume En = 0, 
as well as $ = at the stellar surface, and $ = along 
the last open field line. In addition, the condition E\\ = 
is imposed at an upper boundary, ie, PFF. The OG mod- 
els, which elongate along open field lines, usually have 
upper and lower boundaries along the field lines; an inner 
boundary, which is on the star side; and an outer bound- 
ary, which is located near the light cylinder. In particular, 
the lower boundary is supposed to be the last open field 
line. A standard boundary condition is such that $ = 
on the upper and lower boundaries, $ = and En = 
on the inner boundary, and En = on the outer bound- 
ary. The boundary conditions are very important because 
they affect the structure of the accelerating electric field. 
However, the screening process on each boundary has not 
been well modeled or simulated. Therefore, it is not clear 
whether <f> = along the last open field line or not. 

There is also a problem with the directions of the cur- 
rent flow. The directions of the current along the mag- 
netic field line above the last open field line are different 
in these acceleration models. When \i ■ ft > 0, where /j, 
is the magnetic moment, the direction of the current in 
the PC and SG is inward, and while the current in the 
OG is outward. When /j, ■ ft < 0, the directions are re- 
versed. Looking at various pulse shapes from radio to 
gamma-rays, people sometimes speculate all these types 
of acceleration co-exist. If all of these acceleration models 
exist above the same last open field line, the directions of 
the currents in these models are opposite. Therefore, it 
seems to be impossible that polar-slot gap and OG exist 
above the same last open field lines consistently from the 
theoretical prospective. 

On the other hand, most of the observed radio and 
gamma-ray features of both normal and millisecond pul- 
sars are well reproduced by radio emission from upper 
part of the PC region and gamma-ray emission from the 
SG or OG (eg, Watters et al. 2009; Venter et al. 2009). It 
is required on the observational side that the PC and SG 
and/or OG co-exist. 

To understand the formation mechanism of the PC, SG 
and OG in view of electrodynamics, particle simulation for 
the global magnetosphere would be a powerful way. By 
using GRAPE (Makino et al. 2003), a massively-parallel 
special-purpose computer for astronomical N-body simu- 
lations at National Astronomical Observatory of Japan, 
we performed a simulation particularly aimed at under- 
standing the boundary conditions for the PC, SG and OG; 



the relationship among the models; and the current direc- 
tions. 

The organization of this paper is as follows: Section 2 
describes the method of our simulation; Section 3 shows 
the simulation results and describes the structure of the 
global pulsar magnetosphere; Section 4 discusses the posi- 
tion relationships of the PC, SG, and OG, and summarizes 
the results of our simulation. Our main conclusion, that 
both polar-slot gap and OG can exist in a magnetosphere, 
is also presented in Section 4. 

2. Method of Simulation 

It is necessary to handle formation of the field-aligned 
electric field E\\, so that the simulation must not be the 
usual MHD simulation. A breakthrough can be found by 
a global particle simulation, which can handle cross-field 
drift motions due to radiation drag and particle inertia 
by solving each particle motion. Wada & Shibata (2007, 
hereafter WS) have succeeded in overcoming this heavy 
computation by using GRAPE. Typical use of GRAPE 
is for stellar dynamics. However, it should be noted 
that GRAPE, in which sign-bit is available, can calcu- 
late Coulomb force as well as gravitational force because 
both of them arc represented by the inverse-square law: 
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where the subscripts i, j represent i/jth particle; and 
N,r, e, q, and E are the total particle number, position 
vector, softening parameter, charge, and electric field, re- 
spectively. 

In our simulation, we assumed that a star is an aligned 
rotator and a perfect conductor. To specify the charge 
sign, we chose aligned (/x • ft > 0), not counter aligned 
(fi ■ ft < 0), ie, the poles are negatively charged, while 
the equator is positively charged. The numerical method 
follows WS. Magnetospheric plasmas are represented as 
super-particles. The electric field on the ith particle is 
calculated by the superposition of the vacuum solution 
E v and the space-charge solution E q , where 
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where e r and eg are unit vectors in spherical coordinates; 
R-i Qsysj and c are the stellar radius, system charge and 
speed of light, respectively. To satisfy the boundary con- 
dition for a perfect conductor, the second and third terms 
in (4) representing the effect of mirror charge exist. 
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As an improvement from WS, we take into account the 
modification of the magnetic field from the dipole field. 
By using GRAPE, the magnetic field is calculated by the 
Biot-Savart law: 

/ R \ 3 1 
B(r t ) = bJ — \ (cos6>e r + -sin6>e e ) (5) 

2-f (|r, -nP + f 2 ) 3 / 2 ' k ' 

where Vj is the velocity of the jth particle. 

We made the following assumptions as WS did: (1) par- 
ticles are emitted freely from the stellar surface; (2) the 
particles are subject to the radiation drag force; and (3) 
pair creation occurs where the field-aligned electric field 
En is stronger than a critical value E cr . Since we are look- 
ing for an axisymmetric steady state, we ignored the effect 
of the time variation of the fields. Therefore, the problem 
is electro-static and magneto-static, and whereby GRAPE 
can be used. Our simulation is as follows: starting from 
vacuum around the star, the surface charge is pulled out 
from the star and, as it spreads throughout the magneto- 
sphere, radiation drag takes place; and pairs are created, 
if £j| > E CI . This simulation proceeds in the following 
steps: 

1. Start the calculation from the vacuum around the 
star. 

2. Replace the surface charges on the star by particles. 

3. Compute the electromagnetic field at the particles' 
positions. 

4. Compute the force exerted on the particles, and ad- 
vance the particle positions and velocities by the 
equations of motion. 

5. Create electron-positron pairs where En > E CI . 

6. Remove particles which fall back to the star or leave 
the simulation box through the outer boundary. 

7. Go back to step (2) unless the steady state is estab- 
lished. 

Our simulation is started from the vacuum state. Then, 
the surface charges are pulled out by the induced electro- 
motive force which is much larger than the work function. 
They are emitted as charged particles above the stellar 
surface. The emission is expected to last unless E\\ ^ 0. 
Therefore, we emitted the surface charge a per unit area 
which is defined as En = 4-7r(<7 — aQ])(n ■ B), where <tgj is 
the surface charge density when the corotational electric 
field exists both inside and outside of the star; n is the 
unit normal vector to the stellar surface; and B is the 
unit vector of the magnetic field on the stellar surface. 
The surface charges are emitted ~ 230 times in a rotation 
period. 

We solve the equation of motion, 

5=g,(£ + ftxB) + F rad , (7) 
dr 

where /3,; is the velocity of the zth particle in units of 
the speed of light, p L = jirriifiiC, and 7.; = (1 — /3 2 ) -1 / 2 ; 
F rad = {2/3){q 2 /R 2 c )jt{p l /\p l \) (if ft ~ 1) is the radiation 
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reaction force, and R c is the curvature radius. This ap- 
proximate expression for the radiation drag force can be 
used when the force due to the external field is larger than 
the radiation drag force. 

For pair creation, grid points r.j _j = (ri,6j) are provided 
on the meridional plane. If _Eji is larger than the critical 
value E CI , a\E\\(ri t j)\/E CI pairs whose initial velocities are 
zero are generated in the grid spaces. Here the multiplicity 
a is set to a = 2 and E CI = 0.25£?l, where Bl = n/R^ is the 
light cylinder magnetic field, in the present simulation. 

The volume of the grid spaces is in proportion to r~ 2 . 

so that pair creation rate per unit volume is in propor- 

2 

tion to r . This reproduces the fact that the density of 
target photons supplied from the stellar surface decreases 
in proportion to r 2 . A time interval for the pair creation 
is also the parameter of the calculation. In the present 
calculation, it is such that pairs are created ~ 75 times in 
a rotation period. 

The outer boundary is located at 10 Rj,. 
Magnetospheric plasmas are expressed as super-particles, 
and therefore the exaggerated electric force causes a 
bound motion like positronium. We remove such bound 
pairs because they behave as if they were a neutral 
particle, and they would not affect the structure of the 
magnctosphcrc. In the steady state, creation and loss 
of the particles in the simulation balance with each 
other. Thus, the total charge of the system is determined 
automatically. 

As for normalization, the position vector, magnetic 
field, electric field, and velocity are normalized as follows: 
f = r/R, B = B/B*, E = E/B+, and j3 = v/c, where B* 
is the magnetic field strength on the poles; the variables 
with bar are non-dimensional in this paper. In our simu- 
lation, the strength of the magnetic field is normalized by 
-B*, so that our simulation can treat pulsars which have 
various magnetic field strengths basically. The other ma- 
jor simulation parameters are £7 = R/Ri, = 0.2, fri = 10~ 10 , 
and q = 10 -5 , where Ri, is the light cylinder radius; fh and 
q are the mass and charge of a super-particle, respectively, 
in non-dimensional form. 

The simulation parameters should be chosen so that 
what takes place in the simulation reproduces the actual 
phenomena in the pulsar magnetosphere. At the same 
time, we have a restriction of the computational power. 
Our parameters include the mass and charge of the super- 
particles fh and q, time steps of integration, angular veloc- 
ity Q, frequency of pair creation, critical electric field for 
pair creation E cr , pair multiplicity a, frequency of injec- 
tion from the stellar surface, and softening parameter of 
Coulomb force. It is notable that the total number of par- 
ticles in the simulation box is determined sclf-consistently 
but can be controlled by q and the parameters of pair 
creation. We do not have any grid points for the electro- 
magnetic field. We give typical scales of time and length 
for the simulation, normal pulsars, and MSPs in Table 1 
and 2, where 71 and 72 are the Lorentz factors in units of 
10 and 10 2 , respectively; P-3 and Pi are the stellar rota- 
tion periods in units of 10 -3 sec and 1 sec, respectively; 
B 8 and B i2 are the magnetic field strength in units of 10 8 
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G and 10 12 G, respectively; R§ is the stellar radius in units 
of 10 6 cm; and M% is the multiplicity which is defined as 
M = n/nQj in units of 10 2 . The time step was taken to be 
much smaller than the Larmor period to reproduce drift 
motions due to the electric field, the gradient of the mag- 
netic field, the centrifugal force, and the radiation drag, 
and so on as seen in Table 1 . To consider the screening ef- 
fect and the formation of the field-aligned electric field due 
to inertia, the time step also has to be much smaller than 
the the period of plasma oscillation. This requirement is 
also satisfied as seen in Table 1. According to Table 2, the 
Larmor radii and inertial lengths are enough smaller than 
the typical scale of the system. The ratio of the plasma 
frequency lo p to the gyro frequency ui g changes in such a 
way that oj p /oj g <C 1 at the magnetic poles, io p /uj g < 1 at 
the light cylinder, and u> p /uj g > 1 at the outer boundary 
in both simulation and actual pulsars. Thus, the magne- 
tization property of plasma is also simulated well in the 
simulation and actual pulsars. 

3. Results 

Magnetospheric plasmas are represented by ~ 1 x 10 5 
super-particles. After integration of ~ 10 periods of time, 
the amount of particles produced at the stellar surface 
and the gaps, and the amount of removed particles are 
balanced, so that the total particle number becomes con- 
stant in the simulation box. The net charge of the system 
also converges to a constant value, and therefore we can 
consider the magnetosphere is in a steady state. The top 
left panel of figure 1 shows particle distribution and veloc- 
ity field on the meridional plane in the steady state. Most 
of the positive and negative particles are located around 
the equatorial plane and polar regions, respectively. In ad- 
dition, a near vacuum region appears around the middle- 
altitudes. 

The bottom left panel of figure 1 shows the strength of 
the field-aligned electric field E\\. The strong E\\, which 
is > 0.25£?l, exists around the null charge surface and 
above the poles. In total ~ 3.6 x 10 4 electron-positron 
pairs (super-particle pairs) are generated per rotation. 

The former region is identified as the OG (see WS). 
However, note that the strong E\\ above the poles still has 
artificial effects because it was found that it decreases with 
increasing the frequency of the particle emission on the 
stellar surface. When we changed this emitting frequency 
from ~ 45 times to ~ 230 times per rotation, the 
above the poles decreased from ~ 0.036 to ~ 0.016. This 
indicates that if the emission frequency becomes higher, 
the value of E^ would approach £ji ~ 0. However, it is 
difficult to remove such an artificial effect in the global 
particle simulation due to a limited computational power. 
It is not clear at the moment whether the above the 
poles includes the essential electric field or not. 

The pairs generated in the OG are accelerated in oppo- 
site directions by the E\ | . The positive particles are accel- 
erated outward, and flowed out as a pulsar wind. On the 
other hand, the negative particles are accelerated toward 
the star by the outward-directed E\\. If the star solely 



absorbs the negative particles continuously, it would be 
charged up negatively. To prevent such charging up, the 
negative particles are re-emittcd from the stellar surface 
around the magnetic poles and spread into the magneto- 
sphere. In this way, a current passing through the PC 
originates the negative particles generated in the OG. 

In our simulation, particles are traced up to the distance 
of 10i?L- Most of them flow out beyond the outer bound- 
ary as a pulsar wind, but some negative particles emitted 
from the PC region return to the star. This circulating 
flow can cross the magnetic field due to the -F ra d x B drift. 
In order to understand how modification of the magnetic 
field lines affects the structure of the magnetosphere, we 
have also performed the simulation with the magnetic field 
fixed to the dipole field. In the case of the dipolc field, 
the magnetic flux crossing the equatorial plane beyond i?L 
and 10i?L arc = 0.1 and ^iol = 0.01, respectively. In 
contrast, in the case of the modified magnetic field, the 
values are ~ 0.14 and ^iol ~ 0.03, respectively. The 
closed magnetic fields are pushed out and the magnetic 
flux passing through the outer boundary is increased. This 
indicates that the open flux is increased. The top right 
panel of figure 1 shows the particle distribution and ve- 
locity field in the case of the dipolc field. In comparison 
with the dipole case, it was found that the circulating flow 
diminishes in the case of the modified magnetic field. As 
for the total number of particles going through the outer 
boundary, it is ~ 3.3 x 10 4 per rotation, which corresponds 
to ~ 0.6-ZVgj, where N GJ = (Q 2 R 3 B*)/(2ec) is a current 
of primary particles flowing out from the single polar cap 
in Goldreich & Julian model, in the case of the modified 
magnetic field. 

The non-corotational potential $ represents deviation 
of the electric field from the corotational one. We assume 
$ = on the stellar surface. The bottom right panel of 
figure 1 gives the map of 4>, where if the region is white, 
then $ ps (i.e. the corotation region). In our results, the 
polar regions are reddish due to the potential drop just 
above the stellar surface, while the low-latitude region, 
except for the equatorial closed magnetic flux, is bluish 
due to the outer gap. 

Note that the regions other than the white region in 
the bottom right panel of figure 1 do not always mean 
i?H ^ 0. Even if $ 7^ 0, the field-aligned electric field can 
be zero when B ■ V<1> = 0. On the other hand, the white 
region where 4> = is connected with the star through the 
magnetic field lines on which E\\ =0, and corotates with 
the star. 

There are two distinguishable corotation white regions. 
One is the traditional dead zone with closed magnetic flux 
(the reason why it is rather small shall be discussed later). 
The other is in the middle latitudes extending along the 
field lines. This elongated corotation region separates the 
OG regions with out-going current from the PC regions 
with in-going current. In other words, the middle-latitude 
corotation region is formed where the directions of the cur- 
rents change. Let us call this region "the current-neutral 
zone" . In our numerical experiment, there is no significant 
flow of both positive and negative particles along the mag- 
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Fig. 1. Top left panel: Particle distribution (dots) and velocity fields (arrows) on the meridional plane in the case of the modified 
magnetic field. Red dots and arrows are for positive particles, and blue ones are for negative particles. The star is at the bottom 
left. The dotted lines show the dipole magnetic field lines, and the solid lines show modified magnetic field lines. The geometrical 
scale length is normalized by the stellar radius, and the light cylinder is located at the lateral distance of 5. Top right panel: Particle 
distribution in the case of the dipole magnetic field. Bottom left panel: Color contour of the field-aligned electric field E| | . The 
black line represents where the value of the electric field is equal to E CT . Bottom right panel: Color contour of the non-corotational 
potential <!>. 
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Table 1. Typical periods for the simulation, normal pulsars, and MSPs. 
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6.3 x irr 4 7i 


3.6 x IQ-^^B^ 1 






3.6 x 10- 


16 72B f 


-l 
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5.0 x 10" 1 


2.0 x W- 5 P?B~ 2 1/2 R 6 ' i/2 M 2 


1/2 


2.0 


x 10- 9 P 2 3 B g 
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1.6 x 10- a 7i 


7.8 x 10- 6 72 PfB^Re 3 






7.8 x 10- n 72 




-1 p-3 


Plasma Period at OB 


5.0 x 10 


2.0 x 10- 4 P 2 B 12 1/2 i?e :i/2 M 2 


1/2 


2.0 


x 10- 8 P 2 3 B 8 


1/2 R- 


3/2 M -l/2 


Larmor Period at OB 


1.6 x 10 2 7i 


7.8 x IQ-^PfB^Rf 
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Table 2. Typical lengths for the simulation, normal pulsars, and MSPs. 



Simulation Normal Pulsars (cm) MSPs (cm) 

Stellar Radius 1 10 6 10 6 

Light Cylinder Radius 5 4.8 x 10 9 Pi 4.8 x 10 6 P_ 3 

Larmor Radius at Pole 10 _4 7i 1.7 x lO-^-Bfa 1 1.7 x lO^^B^ 1 

Inertial Length at Pole 5 x 10~ 3 2.0 x W^P^ 2 B~ 2 1/2 M 2 1/2 6.3 x lO^P^ 2 B~ 1/2 M 2 1/2 

Larmor Radius at LC 2.5 x 10~ 2 7i 3.7 x W 4 -y 2 P? B^ 1 R~ 3/2 3.7 x IQ-^P^B^Rq 3/2 

Inertial Length at LC 7.9 x 10~ 2 9.4 x l^PfB^' 2 R~ 3/2 M 2 1/2 9.4P 2 3 B^ 1/2 R~ 3/2 M 2 1/2 

Larmor Radius at OB 2.5 x IO71 3.7 x lO^Pf-B^e 3/2 3.7 x lO^P^B^Pg 3/2 

Inertial Length at OB 7.9 x 10~ 3 9.4 x 10 3 Pf g^P' 372 ]!/' 172 9.4 x lO^P^g' 172 ^ 372 ]!/' 172 



netic field lines in the current-neutral zone. The cloud of 
negative particles is naturally formed there. The current- 
neutral zone is formed in both the dipole and modified 
magnetic field cases, although the locations were slightly 
different because it is formed along the magnetic field 
lines. 

The traditional dead zone is defined by the magnetic 
field lines closed within the light cylinder, but the one in 
our simulation does not reach the light cylinder, and is ob- 
viously small. Additionally, the shape of the disk does not 
follow the magnetic field lines. As shown in figure 2, the 
inner part of the disk corotates with the star, but on the 
other hand the super-rotation, which is caused by the E\\ 
in the gap, took place at the outer part of the disk. The 
Lorentz factor increases to saturate by radiation drag, and 
it is ~ 10 for the super-particle. This Lorentz factor would 
correspond in the actual pulsar to ~ 10 7 P 1 1 ^B X ^R'J A ~~, 
at which the Lorentz force and the radiation drag are equal 
to each other near the light cylinder. For this reason, 
-F ra d x B drift causes cross-field motion within the light 
cylinder, so that the particles leak from the outer part of 
the disk toward the light cylinder. The primary force on 
the particles are the Lorentz force, and therefore shows 
gyro-motion, thus the approximation formula for the ra- 
diation drag force can be applied. The secondary force 
is radiation drag and causes the drift motion. Thus, the 
size of the corotational disk becomes smaller. As for the 
ratio of the magnitude of the electric field to the magni- 
tude of the magnetic field, the electric field is dominant in 
the wedge-shaped region beyond the light cylinder (figure 




b Q.i oa aj a> an a* 0.7 ox o» 1 



Fig. 2. The color-coded azimuthal velocity of the particles 
(circles). The meaning of the colors is shown in the bar just 
below the equatorial plane, where each color indicates the 
corotational speed at each axial distance. 
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Fig. 3. Color contour of the ratio of the magnitude of the 
electric field to the magnitude of the magnetic field, |£!|/|B|. 




Fig. 4. Geometry of the pulsar magnetosphere on the merid- 
ional plane. 

3). This feature, whose apex angle is « 28°, is similar to 
the feature found in the force-free solution by Uzdcnsky 
(2003), whose apex angle is rj 62°. 

4. Discussion and Conclusion 

Let us discuss the structure of the global magnetosphere 
on the basis of our results. In figure 4, we summarize 
electromagnetic properties of each magnetosphcric region. 
There is a possibility that multiple acceleration regions ex- 
ist in a different way from the traditional picture because 
of the current-neutral zone. 



First, we discuss Region A and Region B. In PC and SG 
models, the side boundary condition is that $ = along 
the last open field lines. In our simulation, the side bound- 
ary of the polar flows correspond to the rim of Region B. 
Therefore, the surface of Region B plays the role of the 
side boundary of PC and SG models, instead of the tra- 
ditional dead zone. The boundary condition on the neu- 
tron star surface depends on conditions such as the surface 
temperature, the surface magnetic field strength, and the 
atomic number of matter in the surface layer. In the case 
of most pulsars, one may assume E\\=0 and $ = for the 
boundary conditions on the stellar surface. As described 
in the previous section, the strong E\\ just above the mag- 
netic poles (figure 4) includes the artificial effects, which 
are caused by the limited emitting frequency of the parti- 
cles from the stellar surface. If the emitting frequency is 
increased, the value of the will likely approach zero. 

In space-charge limited flow models, in which Eu = 
and $ = on the stellar surface, particles are acceler- 
ated and emit gamma-rays. Then magnetic pair pro- 
duction occurs, and the PFF is formed. We applied a 
simple calculation method for generating pairs in order 
to avoid heavy computations. However, if we could ap- 
ply a more accurate and realistic calculation, the PFF 
would be properly treated, and the PC and SG would 
appear in a region bounded by the stellar surface and 
Region B. Additionally, we also performed a simulation 
in which pair creation was suppressed above the magnetic 
poles. Even in this simulation, we obtain an steady state 
where the OG is formed. This simulation corresponds to 
a pair-starved polar cap (PSPC) (Muslimov & Harding 
2004; Muslimov & Harding 2009), in which the amount 
of pair creation is not enough to completely screen the 
accelerating electric field. It was found in this simulation 
that the structure of Region B was almost unchanged as 
compared with the previous simulation. This indicates 
that the current-neutral zone is concrete. In summary, 
the acceleration, like polar-slot gap or PSPC, may exist 
in a region bounded by Region B. It is noted that the pair 
creation in the OG is the key to make the magnetosphere 
active even in the aligned rotator, although the diocotron 
instability is also proposed to play a similar role (see, eg, 
Spitkovsky & Arons 2002; Petri 2009). If pair creation 
does not occur anywhere, the dome-disk magnetosphere, 
such as the magnetosphere in Krause-Polstorff & Michel 
(1985), is formed in our simulation. In contrast, if the pair 
creation occurs in the OG, the generated electrons go back 
to the star and are re-emitted around the magnetic poles, 
and then, they flow out to the magnetosphere. In this 
way, pair creation in the OG can change a static magne- 
tosphere into an active magnetosphere without diocotron 
instability. 

Next, let us discuss the region around the OG. The 
corotation Region D, where $ = 0, is identified as a tradi- 
tional dead zone, but Region D is smaller than the tradi- 
tional one and does not reach the light cylinder. Instead, 
Region E exists between Region D and the light cylinder. 
In Region E, E\ =0, but $ ^ as seen in figure 4 because 
the magnetic field lines which pass through Region E go 
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Fig. 5. The diagram of the super- rotation region and outer 
gap. The electric potential charges along on the field line 
due to the outer gap. The potential difference between 
NPin and N&l becomes larger than the corotational one, and 
therefore a stronger Ej_. Thus super-rotation appears. 

through the OG, where Eu ^ 0. The super-rotation arises 
in Region E due to the gap. Note that the super-rotation 
appears even if E\ | = in the disk. This means that rota- 
tion speeds of particles approach the speed of light within 
the light cylinder (see figure 2), and the particles in the 
region can leak out due to the cross-field motion. The 
E\\ in the OG is > 0.25£>l, and the potential drop in the 
OG is ~ 20% of the effective electromotive force in the 
simulation. In the case of actual young pulsars, the pair 
creation rate is higher, so that the strength of the Eu in 
the gap should be smaller. In addition, super-particles 
can cross the field lines easier due to their large masses. 
Therefore, it is thought that the size of Region E in the 
actual magnetospherc is much smaller than that of the 
simulation. 

The azimuthal velocity in the closed field region origi- 
nates from trans-field potential difference via E± x B drift 
motion, which may be given by 

y = cn7— , (8) 

where (f> is the scalar potential. Let us define the coro- 
tating dead zone by \1> < ^/; n , and = (^l/c)^- m on its 
surface (see figure 5). Outside of this corotating dead 
zone, but within the field line tangent to the light cylin- 
der, denoted by vf^, we have a super-rotating region with 
closed field lines (vf^n < * < ^l)- The electric potential 
on the field line "Jl is </>l = (^/c^l on the stellar surface, 
but owing to the potential drop A$ along the field line in 
the gap, the potential on the same field line and on the 
equatorial plane is reduced to 0l — A$. Thus, the poten- 
tial difference between ^> ln and is increased to provide 
super-rotation: 



V<f, 



cw - 



A$) 



7 n i 



> wQ, 



(9) 
(10) 



which indicates that somewhere within the light cylinder, 
the rotation velocity formally exceeds the light speed, say 
at w = vd\. In reality, radiation drag in the azimuthal 
direction drives trans-field motion in a radial direction, 
where just within w\, and beyond w\. In our simulation, 
the thickness of the super-rotation region is artificially en- 
larged because of the somewhat large value of E cr . The 
actual thickness of this region Aw = i?L — tfi , can be es- 
timated by (10). Providing Aw ~ i?L(*in - < 
i? L , we have Aw/R h ~ (A^/^l) 1 / 2 m (/z/10 30 Gcm 3 ) 1 /2p i 
where A$ ~ 10 13 Volt is assumed. In the Crab pulsar, 
Aw ~ 0.03i?L- However, in old pulsars, where $ becomes 
comparable to </>l, the super-rotation region must be sig- 
nificant. 

There would be no pair creation on the magnetic field 
flux between w\ and w- m , and the inner part of the OG. If 
pair creation occurred there, the pairs would not be able 
to cross the closed magnetic field lines and accumulate. If 
so, it is expected that the inner part of the OG is filled 
with the pairs, and the gap would disappear. However, 
this is inconsistent with the initial assumption that the 
gap exists. Thus, the pair-starved region must exist in 
the inner part of the OG. 

In the super-rotational region of the actual pulsar, some 
positrons which generated in the OG flow out at an 
amount approximate to the amount of the GJ current to 
the outside of the magnetosphere by the F x B drift. On 
the other hand, some electrons return to the star, and the 
amount of this current is thought to be that of the GJ 
current, at most. However, the number of pairs generated 
in the OG is ~ 100 times larger than that of the primary 
particles. This indicates that most of the pairs do not be- 
come a current. They might leak as a plasmoid when the 
particle energy density becomes larger than that of the 
magnetic field (see figure 5). We later suggest magnetic 
reconnection at the top of Region E. 

Region E mostly consists of positrons. The positron 
flow is divided into two kinds of flows in the simula- 
tion with the magnetic field fixed to be dipolc. One 
is a flow along the equatorial plane to escape from the 
magnetosphere, and the other one is a flow circulating 
the magnetosphere to return back to the star from the 
poles. However, in the case of the modified magnetic 
field, the circulating component disappears, and most of 
the positrons flow out as a pulsar wind. In the present 
simulation, the mean free path for pair creation is set 
at zero for simplification, and the generated pairs do not 
have initial outward velocities. This approximation works 
well at acceleration regions where pair creation occur ac- 
tively, but does not work well around the acceleration 
regions whose structure would not be reproduced cor- 
rectly in our simulation. Generated pairs are immediately 
charge-separated by the accelerating electric field in the 
gap, and only positrons can escape out of the magneto- 
sphere. In the actual pulsar, the mean free path has a 
finite value, so that gamma-rays can cross the field and 
produce pairs outside the OG. Thus, not only positrons, 
but also electron-positron pairs can escape outside the 
OG. As pointed out by Cheng et al. (1986), gamma-rays 



No. ] 

from the gap always propagate in the convex sides of the 
magnetic field lines, and subsequently, copious electron- 
positron pairs must quench the hcld-aligned electric field 
above the OG (Region C and Region F). In addition to the 
trans-field effect, some pairs generated in the outer part of 
the OG can escape to Region E and Region F due to the 
initial outward velocities. In this way, Region C, Region 
E, and Region F should be filled with quasi-neutral plas- 
mas, so that the condition E\\=Q would be satisfied there. 
If the pair creation rate is made higher than the present 
value so that the current reaches the GJ value, some mag- 
netic field lines would be opened. In that case, the mag- 
netic neutral sheet around the equatorial plane (Region 
G) would be formed. The inner edge of the neutral sheet, 
which is the open-close boundary of the magnetic field 
structure, is called Y-point. Umizaki & Shibata (2010) 
studied the Y-point via an axisymmetric particlc-in-cell 
simulation. They demonstrated that the magnetic recon- 
nection occurred quasi-periodically and thereby plasma 
was heated and accelerated. They also pointed out that 
the electric- field dominant region (\E\ > \B\) was formed 
around the equatorial plane. Our simulation treats steady 
state, so that time-dependent phenomena (eg, the mag- 
netic reconnection) do not appear. However, Region G 
also holds the condition \E\ > \B\ in our simulation (fig- 
ure 3). 

As for the directions of the currents, the direction in 
PC, SG, and PSPC models is inward and the direction 
in the OG model is outward. The mutual location of the 
accelerating regions with opposite currents and how the 
return current is maintained are troublesome problems. 
One possible picture may be indicated in our simulation. 
Note that the inward and outward currents in our simu- 
lation exist along the different field lines separated by the 
current-neutral zone (Region B). 

The size of the OG is determined by balance of the 
amount of supply and loss of pairs. We find in our simula- 
tion that the size of the OG decreases with increasing sup- 
ply of pairs. In the actual pulsars, the supply is influenced 
by the pair creation rate, i.e., the gamma-ray emissiv- 
ity and the mean free path for photon-photon collision to 
make pairs. On the other hand, the loss is due to escaping 
of pairs from the closed field region. However, the loss pro- 
cess is somewhat uncertain, depending on the structure of 
the closed-open boundary and on the super-rotation. The 
super-rotation is affected by the mass-charge ratio of the 
super-particle. As we discussed previously, the actual size 
of the super-rotation region and the gap size in the simu- 
lation can be artificially enlarged. Since the pair creation 
rate is treated as a model parameter and the particle loss 
is merely included quantitatively, it is at the moment dif- 
ficult to discuss actual gap size as function of the pulsar 
parameters such as the period, field strength and X-ray 
flux from the surface. 

Some observational results support the aforementioned 
picture. According to Abdo et al. (2010a), Fermi-LAT 
have detected pulsations from PSR B1509-58 up to 1 GcV 
with a light curve presenting two gamma-ray peaks PI 
and P2, which was observed at phases 0.96 ± 0.01 and 
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0.33 ±0.02, respectively. P2 was consistent with the outer 
magnctosphcre geometry although a sharp cutoff was not 
well explained. Considering the extension to 1 GeV, and 
the magnetic pair and photon-splitting attenuation lim- 
its, PI, which precedes the radio peak, must also origi- 
nate in the outer magnetosphere. However, its phase loca- 
tion could be explained by neither two-pole caustic (TPC) 
models (eg, Dyks & Rudak 2003), which might be real- 
ized in SG acceleration models, nor OG models. PSPC 
models can produce the gamma-ray peak preceding the 
radio peak (Venter et al. 2009), so that PSPC emission 
might be able to explain the phase location of PI. Thus 
to explain both peaks, both pair-starved and non-pair- 
starved (gap) models are required to co-exist (Abdo et al. 
2010a). As mentioned, we also performed the simulation 
without pair creation above the PC, and found that the 
structure around the OG was not affected by eliminating 
pairs in the PCs. This indicates that both pair-starved 
and non-pair-starved gaps exist simultaneously. 

Abdo et al. (2010b) discovered that the distribution of 
the separation of the two gamma-ray peaks appears to 
be bimodal with no strong dependence on pulsar E (or 
age). The A distributions of the separation for each PC, 
TPC, and OG models are calculated by Watters et al. 
(2009). However, the distribution does not show the bi- 
modal feature, but if two acceleration regions exist in a 
magnetosphere, this bimodal feature could be explained. 
The bimodal distribution might indicate that we look at 
different gamma-ray emitting regions from various viewing 
angles. From this viewpoint, we can naturally understand 
the fact that there is no relationship between A and E. 

It is considered that the radio emission from most pul- 
sars come from above the magnetic poles. The radio emis- 
sion is thought to be related to pair creation although the 
mechanism is still not fully understood. Therefore, an ac- 
celeration region to produce pairs is required above the 
magnetic poles. The possibility that gamma-rays orig- 
inate from the PC was ruled out by recent observations 
(eg, Aliu et al. 2008; Abdo et al. 2010b). Nevertheless, ac- 
celeration regions above the magnetic poles, and polar cap 
cascade are required to explain the radio emission. Peak 
separations between radio peaks and gamma-ray peaks 
indicate that acceleration regions are located at different 
regions. These facts also support the idea that at least 
two kinds of acceleration regions exist in a pulsar mag- 
netosphere. Explaining both radio and gamma-ray light 
curves would require both low- and high- altitude accel- 
eration regions. 

In conclusion, the above observational facts strongly 
suggest that both OG and polar-slot gap/PSPC exist. 
From our simulation, these acceleration regions are lo- 
cated at above and below the current-neutral zone. Our 
present simulation is not enough to resolve the structure 
of the PC and SG because the sizes of these regions are 
insignificant compared to the global magnetosphere. As 
a future work, high-accuracy simulations, which have a 
larger number of particles, or follow more realistic pair 
creation processes, are required to discuss the PC and SG 
in more detail. Also, it would be challenging to study the 
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magnetospheric structure in the case of oblique rotators. 
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